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ABSTRACT 



Context. We study the upward propagation of a localized velocity pulse that is initially launched below the transition 

region within the solar atmosphere. The pulse quickly steepens into a shock, which may lead to the formation of 

spicules. 

Aims. We aim to explore the spicule formation scenario in the framework of the rebound shock model. 

Methods. We solve two-dimensional time-dependent magnetohydrodynamic equations numerically to find spatial and 

temporal dynamics of spicules. 

Results. The numerical simulations show that the strong initial pulse may lead to the quasi periodic rising of chro- 

mospheric material into the lower corona in the form of spicules. The periodicity results from the nonlinear wake 

that is formed behind the pulse in the stratified atmosphere. The superposition of raising and falling off plasma 

portions resembles the time sequence of single and double (sometimes even triple) spicules, which is consistent with 

observational findings. 

Conclusions. The two-dimensional rebound shock model may explain the observed speed, width, and heights of type 

I spicules, as well as observed multi-structural and bi-directional flows. The model also predicts the appearance of 

spicules with 3 — 5 min period due to the consecutive shocks. 

Key words. Magnetohydrodynamics (MHD) - Instabilities - Sun: atmosphere 



1. Introduction 

Although spicules were discovered almost 130 years ago 
(Secchi I1877p , they still remain one of mysterious phe- 
nomena in the solar atmosphere. Spicules are usually de- 
tected in chromospheric H a , D3, and Ca II H lines as 
thin and elongated structures in the solar limb. Many au- 
thors consider H a dark mottles as a disc counterpart of 
limb spicules (Suematsu et al. 1995J), however this point 
is still under discussion. Excellent reviews of the general 
properties of spicules (and mottles) were presented al- 
most forty years ago by Beckers (J1968) 11972)) . Suematsu 
(|1998p and Sterling (|2000p published more recent reviews 
of observational and theoretical aspects of spicules, re- 
spectively, while Zaqarashvili & Erdelyi (12Q10|) summa- 
rized observed oscillation events. Very recently, Pasachoff 
et al. (J2009P studied high-resolution dynamics of limb 
spicules by Swedish Solar Telescope (SST) with 0.2 arc sec 
resolution. High resolution observations by Solar Optical 
Telescope (SOT) on board of the recently launched Hinode 
show another type of spicules, which have quite different 
properties than classical limb spicules (De Pontieu et al. 
I2007ap . The classical and newly observed spicules are re- 
ferred to as types I and II, respectively. 

Disc counterparts of limb spicules still are not known 
despite of all the discussion over the years (Tsiropoula et 



al. 11994) Sterling l2000p . However, recent high-resolution 
observations on SST and Hinode/ SOT suggest that active- 
region dynamic fibrils (Hansteen et al. 2006, De Pontieu et 
al. I2007bp and quiet Sun mottles (Rouppe van der Voort 
et al. 2007) have properties similar to type I spicules. On 
the other hand, Rouppe van der Voort et al. (|2009p report 
discovering of disc counterparts to type II spicules through 
observations on SST. 

Despite all the mechanisms of spicule formation pro- 
posed in the literature none of them provides convincing 
explanations of all properties. Spicule formation mecha- 
nisms can be formally divided into three different groups: 
pulses, Alfven waves, and p-mode leakage. 

The main idea of spicule formation by impulsively 
launched perturbations is as follows. Velocity or gas pres- 
sure pulses deposited in the lower atmosphere are steep- 
ened into shocks as a result of the rapid decrease in mass 
density with height. These shocks lift up the transition 
region, producing spicules. Hollweg ([1982)) developed the 
shock model with one-dimensional simulations. In his re- 
bound shock model, Hollweg simulated gas pressure to lift 
the transition region to the observed altitudes. Suematsu 
et al. ()1982p did the same for velocity pulses. Later on, 
Sterli ng et al. ([19901 IT993J) , Cheng (fTM2]) . and Heggland 
et al. ( 200 7) studied the shock models, including radiation 
and heat conduction. 
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Hollweg et al. (j!982p showed that Alfven waves may 
be nonlinear ly coupled to fast magnetoacoustic shocks, 
which may lift up the transition region. Later on, Kudoh 
& Shibata (1999) suggested that the random nonlinear 
Alfvenic pulses may reproduce the spicules at the observed 
heights. On the other hand, the damping of high-frequency 
Alfven waves due to ion-neutral collisions may be respon- 
sible for formation of spicules (Haerendel[1992, James et 
al. 120031) . 

De Pontieu et al. (2004) propose that, as a result of 
fall off of acoustic cut-off frequency, p- modes may be chan- 
neled into the solar corona along inclined magnetic field 
lines. The oscillations then may be steepened into shocks 
producing spicules. De Pontieu et al. (2004) argue that 
the observed quasi 5-min period in spicule appearance is 
associated with the periodicity of p-modes. Several pa- 
pers (Hansteen et al. [20061 De Pontieu et al. 12007b)) then 
studied the formation of spicules due to the leakage of 
photospheric convective motions and oscillations into up- 
per atmosphere. By performing numerical simulations of 
MHD equations the authors show that oscillations and 
flows generated by self consistent two-dimensional (2D) 
convective motions drive spicules. 

Numerous other mechanisms of spicule formation were 
suggested, such as resonant buffeting of anchored mag- 
netic flux tubes by solar granules (Roberts [1979). Detailed 
summary of these mechanisms can be found in the re- 
view paper by Sterling (2000)). Horizontal magnetic flux 
emergence has also been found to be one of many possi- 
ble spicule drivers, the others being convective overshoot, 
collapse of granules, p-modes, and reconnection in either 
the photosphere or chromosphere (Martmez-Sykora et al. 
120091) . 

Despite significant achievements in the above- 
mentioned papers, there are still no existing mechanisms 
that explain the observed double structures of spicules 
(Tanaka [19711 Dara et al. [T9981 Suema tsu et al. I2008J) 
and bi-directional flow (Tsiropoula et al. I1994| Tziotziou 
et al. [20031 [2001 Pasachoff et al I2009J) . Suematsu et al. 
(|2008p and Tsiropoula et al. (1994} have proposed a mag- 
netic reconnection as an explanation of, respectively, the 
double structures and bi-directional flow. However, neither 
cogent numerical simulations nor analytical models have 
been presented so far. Therefore, the problem remains still 
open to discussion. 

In this paper we perform 2D numerical simulations of 
magnetohydrodynamic (MHD) equations and show that 
the 2D rebound shock model may explain both the double 
structures and bi-directional flow in spicules. Unlike pre- 
vious 2D numerical simulations (Hansteen et al. 12006) De 
Pontieu et al. 2007b), which considered a periodic initial 
driver, we launch a single initial pulse and follow its tem- 
poral development. Such an approach allows us to trace 
propagation of the pulse and the internal structure of a 
spicule in detail, which explains double structure and bi- 
directional flow. 

This paper is organized as follows. The main properties 
of limb spicules are summarized in Sect. 2. A numerical 
model is presented in Sect. 3, and the corresponding nu- 
merical results are shown in Sect. 4. This paper concludes 
with a discussion and short summary of the main results 
in Sects. 5 & 6, respectively. 



2. Main properties of limb spicules 

A spicule diameter was estimated from ground-based ob- 
servations to lie within the range of 700 — 2500 km more 
than four decades ago by Beckers (|1968p . Spicules seemed 



to be wider in the Ca II H line than in the H a line (Beckers 
I1972p . Pasachoff et al. (2009]) estimated the mean diameter 
of spicules as 660 ± 200 km using high-resolution observa- 
tions of SST in H a , which fits the old observational data. 
The type II spicules have smaller diameters (< 200 km) 
in the Ca II H line (De Pontieu et al. I2007ap . 
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Fig. 1. Equilibrium configuration of the system: (top panel) 
magnetic field (in units of 1.12 x 10 -3 Tesla) lines are shown 
as arrows, and the mass density g e (y) (colour map plot) is ex- 
pressed in units of 10 -15 kgm -3 ; (bottom panel) temperature 
T e (y) (colour map plot) is drawn in units of 1 MK and block 
boundaries are represented by solid lines. 



While observed in the H a line, type I spicules can reach 
up to 4000 — 12000 km in height from the solar limb with 
a mean value of 7200 ± 2000 km (Pasachoff et al. I2009J) . 
On the other hand, the type II spicules dominate in the 
lower atmospheric layers: they are tallest in coronal holes 
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Fig. 2. Temperature (colour map plot) and velocity (arrows) profiles at £ — 50 s, £ = 190 s, £ = 290 s, £ = 360 s, £ = 400 s, 
t — 500 s, £ = 790 s, and £ = 890 s for (xo = 0, yo = 0.5) Mm and A v = 30 km s _1 . Temperature is drawn in units of 1 MK. The 
arrow below each panel represents the length of the velocity vector, expressed in units of 1 Mm s _1 . The colour bar is common 
to all panels. 
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reaching a level of 5000 km, while they are observed to be 
shorter in quiet and active regions of the solar atmosphere 
(De Pontieu et al. l2007a|) . 

Typical electron temperature and electron density in 
spicules are 15000 - 17000 K and 2-10 n -3.5-10 10 cm" 3 , 
respectively, at altitudes of 4000 — 10000 km above the 
solar surface (Beckers 1968). As a result, spicules are much 
cooler and denser than ambient coronal plasma. 

Spicules arise at a level of about 2000 km with an av- 
erage speed of 25 km s _1 , reach maximum heights, and 
then either fade away or descend back to the photosphere 
with the same speed (Beckers 119681 Pasachoff et al. 2009). 
A typical lifetime of type I spicules is 5-15 min, but some 
spicules may live for longer or shorter periods. Pasachoff 
et al. (2009) find that the distribution of a spicule lifetime 
exhibits a peak at ~ 7 min. On the other hand, the type II 
spicules are characterized by a shorter life time, which is 
about 10-150 s, and they reach higher velocities of 50 — 100 
km s _1 (De Pontieu et al. l2007ap . 

Most spicules show a double thread structure during 
their evolution. This double structure was reported for 
the first time by Tanaka (|1974|) and then by Dara et al. 
(|1998p . Recently, high-resolution observations performed 
by Hinode confirmed that most spicules exhibit the dou- 
ble thread structures (Suematsu et al. I20Q8|) . Suematsu 
et al. (2008) find that the separation of some of the 
double thread spicules vary in time, repeating a single- 
thread phase and the double- thread one. The width of 
each thread and a separation distance between them is a 
few tenths of arc sec on average. 

Another important feature in spicule dynamics is the 
bi-directional flow. This phenomenon has been found to be 
typical of mottles (Tsiropoula et al. 11994} Tziotziou et al. 
12003} I2004p , the base portions showing downward motions 
while the tops of the structures have alternating upward 
and downward motions. On the other hand, Pasachoff et 
al. (J2QQ9|) detected regular, oppositely directed motions at 
lower levels in 40 limb spicules through high-resolution ob- 
servations on SST. Therefore, bi-directional motions seem 
to be typical for spicules. 

3. A numerical model 

Our model system is taken to be composed of a 
gravitationally-stratified solar atmosphere in a 2D space 
x-y. We restrict ourselves to the ideal MHD equations: 

§f + v-u?v) = o, (1) 



<9V 



-Vp+-(V xB) xB 



Qg, 



<9B 



(2) 
(3) 



_ = Vx(VxB), V-B = 0, 
|E + V.(pV) = (l- 7 )pV-V, (4) 

P=^QT. (5) 

m 

Here q is mass density, V is flow velocity, B is the magnetic 
field, p is gas pressure, 7 = 5/3 is the adiabatic index, 
g = (0, —g) is gravitational acceleration of its value g = 
2.72 -10 2 m s -2 , T is temperature, m is mean particle mass 
and k B is the Boltzmann's constant. 



3.1. Initial configuration: equilibrium 

We assume that the solar atmosphere is in static equilib- 
rium (V e = 0) with a force- free magnetic field, i.e., 



(V X B e ) X B e = . 



(6) 



As a result, at this equilibrium the pressure gradient is 
balanced by the gravity force 



Vp e + Qeg = . 



(?) 



Here the subscript e corresponds to equilibrium quantities. 
Using the ideal gas law and the ^-component of hydro- 
static pressure balance indicated by Eq. ([7j), we express 
equilibrium gas pressure and mass density as 



Pe(y) =Po exp 



Here 




Qe(y) 



gMy) ' 



A(y) = k B T e (y)/(mg) 



(8) 
(9) 

(10) 



is the pressure scaleheight, and po denotes the gas pressure 
at the reference level that is chosen at y r = 10 Mm. Note 
that y = corresponds to the base of the photosphere. 

We adopt a smoothed step-function profile for plasma 
temperature 



Te(y) = -T, 



1 + dt + (1 - d t )tanh 



y-yt 

2/w 



(11) 



where d t = T c h/T c with T c h denoting the chromospheric 
temperature at its lower part. The symbol T c corresponds 
to the temperature of the solar corona that is separated 
from the chromosphere by the transition region of its 
width y w = 200 km located at y t = 2 Mm above the 
solar surface. We assume T ch = 15 • 10 3 K and T c = 3 MK. 
For the initial magnetic field, we adopt the magnetic 
field model that was originally described by Priest (1982). 
See also Wasiljew & Murawski (2009) for a recent appli- 
cation of the solar atmosphere model in the context of 
coronal loop oscillations. In this model, we assume that 
Eq. (J6j) is satisfied by a current-free magnetic field 



VxB e 







such that 



B e = V x (Az) . 
Here A denotes the magnetic flux function 
A(x, y) = B A B cos (x/A B )exp[-(y - y r )/A B ] • 



(12) 



The equilibrium magnetic field components (i? ex , B ey ) are 
then given by 

£ ex = #o cos (x/A B ) x exp[-(y - y T )/A B } , (13) 

B ey = -B sin (x/A B ) x exp[-(y - y T )/A B ] (14) 

in addition to B ez = 0. Here, Bo is the magnetic field at 
y = ?/ r , and the magnetic scale- height is Ab = 2L/tt. We 
choose L ~ 30 Mm, which corresponds to the size of a su- 
pergranular cell. This magnetic field is predominantly ver- 
tical at supergranular boundaries (x/A B = n7r/2, where 
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n = ±1,±3, ...), while it reveals a horizontal canopy 
structure at supergranular centres (x/Ab = n7r, where 
n = 0, ±2,...). 

Observations show that type I spicules are mainly 
formed at the supergranular boundaries, where the mag- 
netic field is predominantly vertical. A magnetic field is 
more complicated and is concentrated in tubes in quiet 
regions of the solar photosphere, while at higher altitudes, 
the magnetic tubes merge and produce large-scale struc- 
tures (network and canopy). As we are interested in the 
dynamics of top layers of the solar atmosphere implemen- 
tation of the simple configuration of Eqs. (fT3|) & (fT4|) is 
justified. 

Figure [1] (upper panel) illustrates spatial profiles of 
equilibrium mass density g e (y) (colour map) and mag- 
netic field lines, near supergranular boundaries (arrows). 
The coordinate x=0 corresponds to the central line be- 
tween neighbouring supergranular cells. The mass density 
experiences a sudden fall off with height y at y ~ 0.5 
Mm, and the magnetic field diverges with y. The bottom 
panel shows the spatial profile of the equilibrium temper- 
ature T G (y), which rises abruptly to its coronal value at 
the transition region, y = 2 Mm. 

3.2. Perturbations 

Since we are studying the rebound shock model of spicule 
excitation, we initially perturb (at t = 0) the above equi- 
librium by a Gaussian pulse in the vertical component of 
velocity, viz., 



V y (x, y, t = 0) = A v exp 



(x - x ) 2 + (y- 2/0 ) 



21 



(15) 



Here A v is the amplitude of the pulse, (#o>2/o) is its initial 
position and w denotes its width. 

4. Results of numerical simulations 

Equations ([I])-© are solved numerically using the code 
FLASH (Lee et al. 2009). This code implements a second- 
order unsplit Godunov solver (e.g., Murawski [2002 ) with 
various slope limiters and Riemann solvers, as well as 
adaptive mesh refinement (AMR). We use the mono- 
tonized central slope limiter and the Roe Riemann 
solver (e.g., Toro 1999). We set the simulation box as 
(—3, 3) Mm x (0, 20) Mm and impose boundary conditions 
fixed in time for all plasma quantities in the x- and y- 
directions, while all plasma quantities remain invariant 
along the z-direction. In all our studies we use AMR grid 
with a minimum (maximum) level of refinement set to 3 
(6). The refinement strategy is based on controlling nu- 
merical errors in mass density. This results in an excellent 
resolution of steep spatial profiles and greatly reduces nu- 
merical diffusion at these locations. The initial (at t — 
s) AMR system is displayed in Fig. [T] (lower panel). Every 
block consists of 64 identical numerical cells. 

We launch the velocity pulse of Eq. (14) with yo = 0.5 
Mm (i.e. lower chromosphere) and w = 0.3 Mm. The am- 
plitude of the pulse is taken to be strong enough. We sug- 
gest that the pulse originates at the bottom of the pho- 
tosphere. Therefore its amplitude significantly grows with 
height as a result of the rapid decrease in mass density 



in the photosphere (Carlsson & Stein ll99"7|) . The width of 
the photosphere (500 km) implies about 4 pressure scale- 
heights. Therefore, the amplitude of a velocity signal grows 
by ~ 7 times because of the energy flux conservation in the 
case of vertical propagation. For example, the velocity of a 
granular cell of 1.5 — 2 km s _1 would lead to a flow ampli- 
tude of ~ 10 — 15 km s _1 at height of 500 km. A stronger 
initial pulse at the photospheric level would obviously re- 
sult in a larger amplitude flow at the lower chromosphere. 
Roberts (1979) proposed that the resonant buffeting of 
anchored magnetic flux tubes by granular cells may trig- 
ger strong upward flows of 4 — 5 km s _1 at the base of 
the photosphere. These flows would result in magnitudes 
as large as ~ 30 — 35 km s _1 in the lower chromosphere. 
Other types of activity, such as, a magnetic reconnection, 
may trigger even stronger pulses. Therefore, we consider 
A v > c s , where c s is the chromospheric sound speed that 
is of the order of 10 km s _1 . 

Figure [2] displays the spatial profiles of vertical ve- 
locity (arrows) and plasma temperature resulting from 
the initial velocity pulse that was launched at the point 
(xq = 0,?/o — 0.5) Mm, at which the magnetic field is 
essentially vertical. The upper left panel corresponds to 
t = 50 s. Chromospheric plasma lags behind the shock 
front, which rises up to y = 3.4 Mm. The reason for the 
material being lifted up is the rarefaction of the plasma 
behind the shock front, which leads to low pressure there. 
As a result, the pressure gradient force works against grav- 
ity and forces the chromospheric material to penetrate the 
solar corona. The chromospheric plasma reaches the level 
ofy = 6.5 Mm in the next 100 s. The upper central panel is 
drawn for t = 190 s, which clearly resembles a spicule-like 
structure with the chromospheric temperature and mass 
density. Its width and mean rising speed are 600 — 700 
km and 25 km s _1 , respectively. These values are close 
to the corresponding characteristics of type I spicules. At 
t = 190 s, the plasma already flows downwards because it 
is attracted by gravity. 
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Fig. 3. Timesignatures of velocity V y collected at (x = 0, y = 
2.5) Mm for (x = 0, y = 0.5) Mm and A v = 30 km s" 1 . 
Time t and V y are expressed in units of 1 s and 1 Mm s _1 , 
respectively. 



The next snapshot (upper right panel) was taken at 
t = 290 s. The spicule already subsided to y = 5.2 Mm 
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In the linear approximation, the wake oscillates with the 
acoustic cut-off frequency 
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Fig. 4. Acoustic cutoff wave period P ac (in units of 1 s) vs 
height (in units of 1 Mm) plotted from Eq. (16). 



and the plasma flows downward near the spicule cen- 
tre. However, as the secondary shock lifts up the chro- 
mospheric material (see two small peaks on the left and 
right sides of central spicule), there are upward flows at 
the spicule boundaries. The secondary shock seems to be 
wider than the initial one. Therefore it lifts up the chro- 
mospheric material in the regions out of the initial pulse. 
The central part is overcome by the falling material of the 
spicule, therefore it is dominated by downflows. There are 
clear bi-directional flows in the spicule: downward at the 
centre and upward at the boundaries. This is fully con- 
sistent with the observational findings (Tsiropoula et al. 
[T9941 Tziotziou et al. I2QQ31 [20041 Pasachoff et al. 12009]) . 
After this moment, the central part of the spicule contin- 
ues to subside, while the two small peaks from the previous 
snapshot become taller and exhibit a clear double struc- 
ture in the region of 3.5 Mm < y < 4.5 Mm. This scenario 
is consistent with the observational data of Tanaka (|1974|) . 
Dara et al. ([1998)) . and Suematsu et al. ([2008]) . 

During the subsequent time intervals, the two peaks 
rise, merge, and reach a maximum height of y ~ 7 Mm 
as seen at t = 500 s. Afterwards the material again falls 
back at the central part, while the next shock forces the 
plasma to move upwards at the spicule boundaries. We see 
again two peaks at t = 790 s in the region of 3.5 Mm < 
y < 5.5 Mm. These two peaks rise again and merge to 
resemble a single spicule at t = 890 s (lower right panel). 
The spicule reveals a clear recurrence-like phenomena with 
single and double structures and uni-directional and bi- 
directional flows occurring and disappearing repeatedly in 
time until perturbations finally fade away and the whole 
system relaxes to a static state. 

Owing to their large amplitudes, initial pulses steepen 
rapidly into shocks. Figure [3] illustrates the ^-component 
of velocity that is collected in time at the detection point 
(x = 0,2/ = 2.5) Mm for the initial pulse amplitude of 
A v = 30 km s _1 . The arrival of the shock front at the 
detection point, y = 2.5 Mm, is clearly seen at t ~ 30 
s. The second shock front reaches the detection point at 
t ~ 320 s, i.e., after ~ 5 min (which fits the observed 5- 
min periodicity in the spicule appearance). This secondary 
shock results from the nonlinear wake, which lags behind. 



^ 'an — 




(16) 



Figure 2] displays the acoustic cut-off period P ac = 
27r/^ ac of the model atmosphere vs height. In the low 
chromosphere, P ac ~ 200 s and then quickly increases 
towards the corona. It is noteworthy that Eq. ([T6]) is 
obtained from the linear analysis, while the nonlinear 
description may change the oscillation period of wake. 
Figure [5] shows the ^-component of velocity that is col- 
lected in time at the detection point (x = 0, y = 2.5) Mm 
for the weaker initial pulse with its amplitude A v = 5 
km s -1 . Indeed, the consequent shocks arrive with signif- 
icantly shorter time intervals than in the case of A v = 30 
km s _1 . The second shock front arrives in t ~ 200 s af- 
ter the passage of the first one. The interval between next 
shocks increases up to 250 s. Thus, the time interval be- 
tween consecutive shocks depends on the amplitude of ini- 
tial pulse, and it is longer for stronger pulses. 




Fig. 5. Timesignatures of V y collected at (x — 0, y = 2.5) Mm 
for (xq = 0,yo = 0.5) Mm, when A v = 5 km s _1 . Time t and 
V y are expressed in units of 1 s and 1 Mm s _1 , respectively. 



This scenario composes the building block of the ID re- 
bound shock model of Hollweg (|1982p , who proposed that 
the secondary shock (or the rebound shock) lifts the tran- 
sition region higher than the first shock, thereby resulting 
in a spicule appearance at observed heights. The process 
is well studied in the framework of ID numerical simu- 
lations. However, our 2D numerical simulations introduce 
interesting new features in comparison to the ID rebound 
shock scenario. 

In the case of the simulations whose results are pre- 
sented in Fig. [2j the chromospheric plasma rises verti- 
cally, because the pulse was initially launched along ver- 
tical magnetic field lines. However, most spicules are in- 
clined to the vertical. Actually, magnetic field is vertical 
only at the centre of a chromospheric network, while it 
is inclined in its surroundings. To study what happens 
around a chromospheric network, we launched the initial 
pulse along the inclined magnetic field. Figure [6] displays 
the spatial profiles of vertical velocity (arrows) and plasma 
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Fig. 6. Temperature (colour map plot) and velocity (arrows) profiles at t = 20 s, t = 50 s, t = 110 s, t = 180 s, t — 290 s, 
£ = 330 s, t = 360 s, and t = 500 s for (xo = 0, 2/0 = 0.5) Mm and A v = 30 km s _1 . Temperature is drawn in units of 1 MK. The 
arrow below each panel represents the length of the velocity vector, expressed in units of 1 Mm s _1 . The colour bar is common 
to all panels. 
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temperature at consecutive times after the initial velocity 
pulse was launched at the point (xo = 1, 2/o = 0.5) Mm, at 
which the magnetic field is slightly inclined (see top panel 
of Fig. [1]). The left top panel shows the profiles at t = 20 
s: the pulse is still propagating within the chromosphere. 
The next snapshot (central top panel) is taken at t = 50 
s. The shock already propagates within the corona with 
high speed, so it reaches up to the level of y = 6.4 Mm at 
this time. On the other hand, the chromospheric material 
lags behind the shock along the inclined trajectory and 
with a much slower speed reaching the altitude of y ~ 3.2 
Mm. The next two snapshots show how the spicule rises 
to the height of ~ 6.5 Mm during the next 130 s. The 
spicule is inclined to the vertical, and it follows the mag- 
netic field structure. The next two snapshots show how 
the spicule subsides, while the chromospheric plasma rises 
at its boundaries due to the secondary shock. The snap- 
shot at t = 330 s shows a triple structure between y = 3.5 
Mm and y = 4 Mm. The last two snapshots show how the 
two lateral peaks rise further with height. First, they form 
the double structure at t = 360 s in the region of 3.5 Mm 
< y < 4.5 Mm. Subsequently, the peaks merge and form a 
single inclined structure, which reaches the level of y = 7 
Mm. 

Figure [7] illustrates the ^/-component of velocity that is 
collected in time at the detection point (x = 1.15, y = 2.5) 
Mm. This figure is very similar to Fig. [3j because both of 
them reveal that the secondary shock arrives after ~ 5 
min. This is because the magnetic field is almost verti- 
cal in the chromosphere near supergranular boundaries, 
where the initial pulses were launched, while it becomes 
more inclined in higher coronal regions (see Fig. [1]). The 
periodicity of the nonlinear wake, which determines the 
time interval between consecutive shocks, is formed in the 
lower region. Therefore, we expect that the period of con- 
secutive shocks does not depend significantly on x in the 
vicinity of a supergranular boundary. However, as we al- 
ready discussed above, the period of the nonlinear wake 
significantly depends on the amplitude of the initial pulse. 
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Fig. 7. Timesignatures of V y collected at (x = 1.15, y = 2.5) 
Mm for (xo = l,j/o = 0.5) Mm and A v = 30 km s" 1 . Time t 
and V y are expressed in units of 1 s and 1 Mm s _1 , respectively. 



5. Discussion 

Our numerical simulations show that the 2D rebound 
shock model leads to the observed dynamics of type 
I spicules. The observed rising speed, width and mean 
length of spicules can be nicely reproduced with a suffi- 
ciently stronger initial velocity pulse. The ID models may 
also produce the same properties. However, the advantage 
of the 2D model is that it may produce other observed 
properties of spicules, such as a double (sometimes multi) 
structure (Tanaka I1974[ Dara et al. [1998, Suematsu et 
al. I2008P and bi-directional flows (Tsiropoula et al. 11994} 
Tziotziou et al. [20031 [2001 Pasachoff et al. [2009]) . This 
is clearly seen in Figs. [2] & [6j The multi-structure and 
bi-directional flows arise when the chromospheric mate- 
rial, which was raised by the initial shock, falls back due 
to the gravity, while the secondary shock is lifting up 
another portion of chromospheric plasma. The chromo- 
spheric plasma is raised up by the pressure gradient behind 
the shock. Therefore, the shock propagates with a higher 
speed that is close to the coronal sound speed, while the 
chromospheric plasma is rising with lower velocity, as ob- 
served in spicules. 

The mean rising speed and maximum height of spicules 
significantly depend on a pulse amplitude A v . Figure [8] il- 
lustrates the dependence of a rising speed (crosses), maxi- 
mum height measured from the transition region (x's), and 
width of spicules (asterisks) on A v . The speeds, height, 
and width are expressed in units of 1 km s _1 , 1 Mm, 
and 0.1 Mm, respectively. The mean speed and maximum 
height strongly depend on A v and both grow with A v . On 
the other hand, the width varies only slightly with A v . It 
turns out that the observed upward speed of a spicule in 
the corona (25 km s _1 ) and its mean height of ~ 7 Mm 
(evaluated from the photospheric level) can be achieved 
for A v = 30 km s _1 . The initial amplitude of velocity 
pulse seems to be quite large for the low chromosphere. 
Pulses from granular motions or p-modes can not reach 
such high amplitudes at ~ 0.5 Mm level; however, res- 
onant buffeting of granules on anchored magnetic tubes 
(Roberts 1979) and/or magnetic reconnection may pro- 
duce such strong pulses. 

The 2D rebound shock model suggests a quasi-periodic 
rise of spicules with the period of ~ 5 min, which is consis- 
tent with observations (Kulidzanishvili & Nikolskv 11978} 
De Pontieu et al. [2001 Xia et al. [2005]) . In our simula- 
tions, this periodicity results from a nonlinear wake that 
is formed behind a leading pulse rather than from p-mode 
leakage. The period of the wake strongly depends on the 
amplitude of the initial pulse and is longer for stronger 
pulses (see Figs. [3]&[5j). The time interval between first 
and secondary shocks is 200 s and 300 s for the pulses 
with initial amplitudes of 5 km s _1 and 30 km s _1 , respec- 
tively. Thus, the periodicity of spicule appearance should 
be 3 — 5 min and even longer, depending on a pulse ampli- 
tude. This is consistent with high-resolution observations 
on SST (Rouppe van der Voort et al. [2007]) . 

In our simulations, the initial pulse was launched 
within the chromosphere, where photospheric magnetic 
tubes already merged and produced a larger scale network 
magnetic field. The field is nearly vertical at supergranular 
boundaries and has an almost horizontal canopy structure 
above the supergranular cells. Therefore, our arcade model 
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Fig. 8. Simulated rising speed (crosses, in units of 1 km s _1 ), 
maximum height (x, in 1 Mm), and width of spicules (asterisks, 
in 100 km) vs initial amplitude of the velocity pulse A v (in 1 
km s _1 ) for (xo — O^yo — 0.5) Mm. Height is measured from 
the transition region, therefore spicule heights from the solar 
surface can be produced by adding of 2 Mm. 



of magnetic field is a good approximation in this region. 
The network structure of the magnetic field is probably the 
reason spicules mostly appear near supergranular bound- 
aries; the horizontal magnetic field prohibits raising chro- 
mospheric material, therefore spicules may only appear at 
the boundaries, where the magnetic field is predominantly 
vertical. It would be also interesting to launch pulses at 
the photospheric level; however, it requires considering a 
more complicated magnetic field structure as it is concen- 
trated in tubes there. This should be the subject of future 
studies. 

We must note here that the model is quite sim- 
plified as the equations do not include such terms as 
plasma partial ionization, thermal conduction, and radi- 
ation. The presence of neutral atoms in plasma is known 
to change its dynamical and physical properties signifi- 
cantly (Braginskii 1965, Haerendel 1992, Khodachenko & 
Zaitsev 2002, Khodachenko et al. 2004). Different interac- 
tion of electrons, ions, and neutral atoms with the mag- 
netic field and with each other causes the main specifics 
of the partially ionized plasma MHD, which differs signif- 
icantly from the fully ionized plasma case. However, the 
effect is more pronounced across the magnetic field, so 
it may have little influence on the rebound shock model. 
On the other hand, radiation and thermal conduction may 
change the scen ario a lot fSterling |T990l[T993l Cheng[1992j 
Heggland et al. I20Q7|) : therefore, they should be included 
in future models in order to have a comprehensive descrip- 
tion of spicule formation process. 

6. Conclusions 

We performed 2D numerical simulations of velocity pulse 
propagation in the solar atmosphere. The strong longi- 
tudinal pulse was launched at the lower part of the chro- 
mosphere, and its consequent propagation along transition 
region into the corona was traced. The pulse quickly steep- 
ens into shock, which influences the dynamics of the up- 



per chromosphere/transition region; namely, the rarefied 
tail behind the shock front leads to the pressure gradient 
above the transition region, which lifts up the chromo- 
spheric material in the form of spicule. The strong initial 
pulse may raise the chromospheric plasma up to observed 
heights (6 — 7 Mm) and with observed speed (25 km s _1 ). 
The spicule begins to falls back after some time, while the 
secondary shock (rebound shock) lifts up another portion 
of the chromospheric material. The superposition of falling 
off and rising plasma portions produces the observed dou- 
ble structure and bi-directional flows in spicules. The sim- 
ulated spicule exhibits a temporal development with dou- 
ble and single structures (sometime even triple) that are 
akin to the observational data. The model predicts quasi- 
periodic raising of spicules with the period of the nonlinear 
wake that is formed behind a leading pulse in the strati- 
fied atmosphere. The periodicity strongly depends on the 
initial amplitude of the pulse and can be in the range of 
3 — 5 min and even longer. We believe that future sophis- 
ticated numerical simulations may give a more complete 
picture of spicule formation in the frame of the 2D and 
3D rebound shock models. 
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